9 Reptiles

9.1 Chalcides striatus

9.1.1 Retrieve data

species="Chalcides striatus"
genus=strsplit(species, " ")[[1]][1] %>% tolower()

sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
    rename(dataset=Dataset) %>%
    filter(Species == species)


read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0122/DMB0122_counts.tsv.gz") %>%
  rename(genome = 1)

genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0122/DMB0122_mag_info.tsv.gz") %>%
  rename(genome = 1, length=mag_size)

genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0122/DMB0122_coverage.tsv.gz") %>%
  rename(genome = 1)

download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0122/DMB0122.tree.gz", "data/DMB0122.tree.gz")
genome_tree <- read_tree("data/DMB0122.tree.gz")

9.1.2 Filter data

#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
  mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
  mutate(across(-1, ~ . * read_counts[[cur_column()]]))

# Transform  to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

9.1.3 Community barplot

# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

# Stacked barplot
genome_counts %>%
    mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE)))  %>% #apply TSS nornalisation
    pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
    left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
    left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
    mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
    ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
          geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
          scale_fill_manual(values=phylum_colors) +
          labs(y = "Relative abundance") +
          facet_nested(. ~ Sample + Extraction, scales="free_x") +
          guides(fill = guide_legend(ncol = 3)) +
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                          axis.title.x = element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
                legend.position="none",
                legend.title=element_blank())

ggsave(paste0("figures/barplot_",genus,".pdf"))

9.1.4 Alpha diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=0) %>%
            t() %>%
            as.data.frame() %>%
            rename(richness=1) %>%
            rownames_to_column(var="dataset")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1) %>%
            t() %>%
            as.data.frame() %>%
            rename(neutral=1) %>%
            rownames_to_column(var="dataset")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1,tree=genome_tree) %>%
            t() %>%
            as.data.frame() %>%
            rename(phylogenetic=1) %>%
            rownames_to_column(var="dataset")

# Merge alpha diversities
alpha_diversity <- richness %>%
      full_join(neutral,by=join_by(dataset==dataset)) %>%
      full_join(phylogenetic,by=join_by(dataset==dataset))

# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))

# Print alpha diversity
alpha_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
  group_by(Extraction) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_lcgg8gw7culp8p96l11v
Extraction richness neutral phylogenetic
DREX 60.5 31.68716 5.314585
EHEX 64.5 34.15145 5.180680
ZYMO 74.0 37.04400 5.254353

9.1.5 Beta diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C", out="pair")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C", out="pair")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")

# Merge beta diversities
beta_diversity <- richness %>%
      full_join(neutral,by=c("first", "second")) %>%
      full_join(phylogenetic,by=c("first", "second")) %>%
      rename(richness=C.x, neutral=C.y, phylogenetic=C)

# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))

# Print beta diversities
beta_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
  mutate(relationship = case_when(
    Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
    Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
    TRUE ~ NA_character_  # Handle other cases if needed
  )) %>%
  filter(relationship != is.na(relationship)) %>%
  group_by(relationship) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_shf7pcijp5l0oh4lk4vg
relationship richness neutral phylogenetic
inter 0.4109306 0.45373698 0.09334837
intra 0.1286148 0.06449944 0.02626687

9.1.6 Variance partitioning

richness <-genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
      write_tsv(paste0("results/var_",genus,".tsv"))

9.2 Natrix astreptophora

9.2.1 Retrieve data

species="Natrix astreptophora"
genus=strsplit(species, " ")[[1]][1] %>% tolower()

sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
    rename(dataset=Dataset) %>%
    filter(Species == species)


read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0125/DMB0125_counts.tsv.gz") %>%
  rename(genome = 1)

genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0125/DMB0125_mag_info.tsv.gz") %>%
  rename(genome = 1, length=mag_size)

genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0125/DMB0125_coverage.tsv.gz") %>%
  rename(genome = 1)

download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0125/DMB0125.tree.gz", "data/DMB0125.tree.gz")
genome_tree <- read_tree("data/DMB0125.tree.gz")

9.2.2 Filter data

#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
  mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
  mutate(across(-1, ~ . * read_counts[[cur_column()]]))

# Transform  to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

9.2.3 Community barplot

# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

# Stacked barplot
genome_counts %>%
    mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE)))  %>% #apply TSS nornalisation
    pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
    left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
    left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
    mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
    ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
          geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
          scale_fill_manual(values=phylum_colors) +
          labs(y = "Relative abundance") +
          facet_nested(. ~ Sample + Extraction, scales="free_x") +
          guides(fill = guide_legend(ncol = 3)) +
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                          axis.title.x = element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
                legend.position="none",
                legend.title=element_blank())

ggsave(paste0("figures/barplot_",genus,".pdf"))

9.2.4 Alpha diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=0) %>%
            t() %>%
            as.data.frame() %>%
            rename(richness=1) %>%
            rownames_to_column(var="dataset")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1) %>%
            t() %>%
            as.data.frame() %>%
            rename(neutral=1) %>%
            rownames_to_column(var="dataset")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1,tree=genome_tree) %>%
            t() %>%
            as.data.frame() %>%
            rename(phylogenetic=1) %>%
            rownames_to_column(var="dataset")

# Merge alpha diversities
alpha_diversity <- richness %>%
      full_join(neutral,by=join_by(dataset==dataset)) %>%
      full_join(phylogenetic,by=join_by(dataset==dataset))

# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))

# Print alpha diversity
alpha_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
  group_by(Extraction) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_5syfync1wl0f50q827kk
Extraction richness neutral phylogenetic
DREX 37.5 16.30449 3.270384
EHEX 29.5 15.34095 3.211162
ZYMO 24.0 13.43764 3.065722

9.2.5 Beta diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C", out="pair")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C", out="pair")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")

# Merge beta diversities
beta_diversity <- richness %>%
      full_join(neutral,by=c("first", "second")) %>%
      full_join(phylogenetic,by=c("first", "second")) %>%
      rename(richness=C.x, neutral=C.y, phylogenetic=C)

# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))

# Print beta diversities
beta_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
  mutate(relationship = case_when(
    Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
    Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
    TRUE ~ NA_character_  # Handle other cases if needed
  )) %>%
  filter(relationship != is.na(relationship)) %>%
  group_by(relationship) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_qkf7b7pmvbetftm642id
relationship richness neutral phylogenetic
inter 0.5386629 0.65025916 0.18128644
intra 0.2756370 0.04665216 0.01166948

9.2.6 Variance partitioning

richness <-genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
      write_tsv(paste0("results/var_",genus,".tsv"))

9.3 Podarcis muralis

9.3.1 Retrieve data

species="Podarcis muralis"
genus=strsplit(species, " ")[[1]][1] %>% tolower()

sample_metadata <- read_tsv(paste0("data/metadata.tsv")) %>%
    rename(dataset=Dataset) %>%
    filter(Species == species)


read_counts <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0128/DMB0128_counts.tsv.gz") %>%
  rename(genome = 1)

genome_metadata <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0128/DMB0128_mag_info.tsv.gz") %>%
  rename(genome = 1, length=mag_size)

genome_coverage <- read_tsv("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0128/DMB0128_coverage.tsv.gz") %>%
  rename(genome = 1)

download.file("https://sid.erda.dk/share_redirect/BaMZodj9sA/DMB/DMB0128/DMB0128.tree.gz", "data/DMB0128.tree.gz")
genome_tree <- read_tree("data/DMB0128.tree.gz")

9.3.2 Filter data

#Filter by coverage
min_coverage=0.3
read_counts_filt <- genome_coverage %>%
  mutate(across(where(is.numeric), ~ ifelse(. > min_coverage, 1, 0))) %>%
  mutate(across(-1, ~ . * read_counts[[cur_column()]]))

# Transform  to genome counts (non-filtered)
readlength=150
genome_counts <- read_counts %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

# Transform to genome counts (coverage-filtered)
readlength=150
genome_counts_filt <- read_counts_filt %>%
  mutate(across(where(is.numeric), ~ . / (genome_metadata$length / readlength) ))

9.3.3 Community barplot

# Retrieve EHI taxonomy colors
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    select(colors) %>%
    pull()

# Stacked barplot
genome_counts %>%
    mutate_at(vars(-genome), funs(./sum(., na.rm = TRUE)))  %>% #apply TSS nornalisation
    pivot_longer(-genome, names_to = "dataset", values_to = "count") %>%
    left_join(., genome_metadata, by = join_by(genome == genome)) %>% #append taxonomy
    left_join(., sample_metadata, by = join_by(dataset == dataset)) %>% #append taxonomy
    mutate(Extraction=factor(Extraction,levels=c("ZYMO","DREX","EHEX"))) %>%
    ggplot(., aes(x=dataset,y=count,fill=phylum, group=phylum))+ #grouping enables keeping the same sorting of taxonomic units
          geom_bar(stat="identity", colour="white", linewidth=0.1)+ #plot stacked bars with white borders
          scale_fill_manual(values=phylum_colors) +
          labs(y = "Relative abundance") +
          facet_nested(. ~ Sample + Extraction, scales="free_x") +
          guides(fill = guide_legend(ncol = 3)) +
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
                          axis.title.x = element_blank(),
                panel.background = element_blank(),
                panel.border = element_blank(),
                panel.grid.major = element_blank(),
                panel.grid.minor = element_blank(),
                axis.line = element_line(linewidth = 0.5, linetype = "solid", colour = "black"),
                legend.position="none",
                legend.title=element_blank())

ggsave(paste0("figures/barplot_",genus,".pdf"))

9.3.4 Alpha diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=0) %>%
            t() %>%
            as.data.frame() %>%
            rename(richness=1) %>%
            rownames_to_column(var="dataset")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1) %>%
            t() %>%
            as.data.frame() %>%
            rename(neutral=1) %>%
            rownames_to_column(var="dataset")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hilldiv(.,q=1,tree=genome_tree) %>%
            t() %>%
            as.data.frame() %>%
            rename(phylogenetic=1) %>%
            rownames_to_column(var="dataset")

# Merge alpha diversities
alpha_diversity <- richness %>%
      full_join(neutral,by=join_by(dataset==dataset)) %>%
      full_join(phylogenetic,by=join_by(dataset==dataset))

# Write alpha diversities
alpha_diversity %>% write_tsv(paste0("results/alpha_",genus,".tsv"))

# Print alpha diversity
alpha_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(dataset==dataset)) %>%
  group_by(Extraction) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_uadht8ha3mkoybb01abw
Extraction richness neutral phylogenetic
DREX 80.5 33.09492 4.689994
EHEX 73.0 32.71816 4.716865
ZYMO 78.5 34.42602 4.745518

9.3.5 Beta diversity

#Calculate Hill numbers
richness <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C", out="pair")

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C", out="pair")

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C", out="pair")

# Merge beta diversities
beta_diversity <- richness %>%
      full_join(neutral,by=c("first", "second")) %>%
      full_join(phylogenetic,by=c("first", "second")) %>%
      rename(richness=C.x, neutral=C.y, phylogenetic=C)

# Write beta diversities
beta_diversity %>% write_tsv(paste0("results/beta_",genus,".tsv"))

# Print beta diversities
beta_diversity %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(first==dataset)) %>%
  left_join(sample_metadata %>% select(dataset,Extraction,Sample),by=join_by(second==dataset)) %>%
  mutate(relationship = case_when(
    Extraction.x == Extraction.y & Sample.x != Sample.y ~ "inter",
    Extraction.x != Extraction.y & Sample.x == Sample.y ~ "intra",
    TRUE ~ NA_character_  # Handle other cases if needed
  )) %>%
  filter(relationship != is.na(relationship)) %>%
  group_by(relationship) %>%
  summarise(richness=mean(richness), neutral=mean(neutral), phylogenetic=mean(phylogenetic)) %>%
  tt()
tinytable_cpqpdqrpg6tkteayq32q
relationship richness neutral phylogenetic
inter 0.2864067 0.41508917 0.122770382
intra 0.0446703 0.00905735 0.003367558

9.3.6 Variance partitioning

richness <-genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=0, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

neutral <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

phylogenetic <- genome_counts_filt %>%
            column_to_rownames(var="genome") %>%
            select(where(~!all(. == 0))) %>%
            hillpair(.,q=1, tree=genome_tree, metric="C") %>%
            as.dist() %>%
            adonis2(formula = . ~ Sample + Extraction, data = sample_metadata %>% arrange(match(dataset,colnames(genome_counts_filt))), permutations = 999, method = "bray", sqrt.dist = TRUE) %>%
            broom::tidy() %>%
            filter(term == "Extraction") %>%
            select(R2) %>% pull()

tibble(species=rep(species,3),metric=c("richness","neutral","phylogenetic"),r2=c(richness,neutral,phylogenetic)) %>%
      write_tsv(paste0("results/var_",genus,".tsv"))